home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / pmlsrc23.zoo / pmlsrc / cexp.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-03-19  |  4.2 KB  |  199 lines
  1. /************************************************************************
  2.  *                                    *
  3.  *                N O T I C E                *
  4.  *                                    *
  5.  *            Copyright Abandoned, 1987, Fred Fish        *
  6.  *                                    *
  7.  *    This previously copyrighted work has been placed into the    *
  8.  *    public domain by the author (Fred Fish) and may be freely used    *
  9.  *    for any purpose, private or commercial.  I would appreciate    *
  10.  *    it, as a courtesy, if this notice is left in all copies and    *
  11.  *    derivative works.  Thank you, and enjoy...            *
  12.  *                                    *
  13.  *    The author makes no warranty of any kind with respect to this    *
  14.  *    product and explicitly disclaims any implied warranties of    *
  15.  *    merchantability or fitness for any particular purpose.        *
  16.  *                                    *
  17.  ************************************************************************
  18.  */
  19.  
  20.  
  21. /*
  22.  *  FUNCTION
  23.  *
  24.  *    cexp   complex double precision exponential
  25.  *
  26.  *  KEY WORDS
  27.  *
  28.  *    cexp
  29.  *    complex functions
  30.  *    machine independent routines
  31.  *    math libraries
  32.  *
  33.  *  DESCRIPTION
  34.  *
  35.  *    Computes double precision complex exponential of
  36.  *    a double precision complex argument.
  37.  *
  38.  *  USAGE
  39.  *
  40.  *    COMPLEX cexp (z)
  41.  *    COMPLEX z;
  42.  *
  43.  *  REFERENCES
  44.  *
  45.  *    Fortran 77 user's guide, Digital Equipment Corp. pp B-13
  46.  *
  47.  *  PROGRAMMER
  48.  *
  49.  *    Fred Fish
  50.  *    Tempe, Az 85281
  51.  *    (602) 966-8871
  52.  *
  53.  *  INTERNALS
  54.  *
  55.  *    Computes complex exponential of z = x + j y from:
  56.  *
  57.  *        1.    r_cexp = exp(x) * cos(y)
  58.  *
  59.  *        2.    i_cexp = exp(x) * sin(y)
  60.  *
  61.  *        3.    cexp(z) = r_cexp + j i_cexp
  62.  *
  63.  */
  64.  
  65. #if !defined (__M68881__) && !defined (sfp004)
  66.  
  67. #include <stdio.h>
  68. #include <math.h>
  69. #include "pml.h"
  70.  
  71.  
  72. COMPLEX cexp (z)
  73. COMPLEX z;
  74. {
  75.     COMPLEX result;
  76.     double expx;
  77.  
  78.     expx = exp (z.real);
  79.     result.real = expx * cos (z.imag);
  80.     result.imag = expx * sin (z.imag);
  81.     return (result);
  82. }
  83.  
  84. #else    !defined (__M68881__) && !defined (sfp004)
  85. #ifdef    ERROR_CHECK
  86. __asm("
  87. .text
  88. .even
  89. _funcname:
  90.     .ascii \"cexp\\0\"
  91.     .even
  92. ");
  93. #endif    ERROR_CHECK
  94. #endif    !defined (__M68881__) && !defined (sfp004)
  95.  
  96. #ifdef    __M68881__
  97.  
  98. __asm("
  99.     .text 
  100.     .even
  101.     .globl _cexp
  102. _cexp:
  103.     fmovex    fp2,sp@-    | 12 Bytes
  104.     movel    a1,d0        | save a1 as return value
  105.     fetoxd    a7@(16),fp0    | exp( z.real )
  106.     fmoved    a7@(24),fp2
  107.     fcosx    fp2,fp1
  108.     fsinx    fp2,fp2
  109.     fmulx    fp0,fp1        |
  110.     fmulx    fp0,fp2        |
  111.     fmoved    fp1,a1@        | fetch result.real
  112.     fmoved    fp2,a1@(8)    | fetch result.imag
  113.     fmovex    sp@+,fp2
  114. ");    /* end asm    */
  115. #endif    __M68881__
  116.  
  117. #ifdef    sfp004
  118. __asm("
  119. | double precision floating point stuff for Atari-gcc using the SFP004
  120. | developed with gas
  121. |
  122. | double precision complex exponential function
  123. |
  124. | M. Ritzert (mjr at dmzrzu71)
  125. |
  126. | 12.10.1990
  127. |
  128. | addresses of the 68881 data port. This choice is fastest when much data is
  129. | transferred between the two processors.
  130.  
  131. comm =     -6
  132. resp =    -16
  133. zahl =      0
  134.  
  135. | waiting loop ...
  136. |
  137. | wait:
  138. | ww:    cmpiw    #0x8900,a1@(resp)
  139. |     beq    ww
  140. | is coded directly by
  141. |    .long    0x0c688900, 0xfff067f8
  142. | and
  143. | www:    tst.b    a1@(resp)
  144. |    bmi.b    www
  145. | is coded by
  146. |    .word    0x4a68,0xfff0,0x6bfa        | test
  147.  
  148.     .text; .even
  149.     .globl _cexp
  150. _cexp:
  151.     movel    a1,d0                | save a1 as return value
  152.     lea    0xfffa50,a0            | fpu address
  153.     movew    #0x5410,a0@(comm)        | exp()
  154.     cmpiw    #0x8900,a0@(resp)        | check
  155.     movel    a7@(4),a0@            | load arg_hi
  156.     movel    a7@(8),a0@            | load arg_low
  157.  
  158. |    movew    #%0101 0101 0011 0001,a0@(comm)    | sincos: sin -> fp2
  159.     movew    #0x5531,a0@(comm)        | sincos: sin -> fp2
  160.     .long    0x0c688900, 0xfff067f8
  161.     movel    a7@(12),a0@            | load arg_hi
  162.     movel    a7@(16),a0@            | load arg_low
  163.  
  164. |    movew    #%0000 0000 1010 0011,a0@(comm)    | mul fp0 -> fp1
  165.     movew    #0x00a3,a0@(comm)        | mul fp0 -> fp1
  166.     .word    0x4a68,0xfff0,0x6bfa        | test
  167.  
  168. |    movew    #%0000 0001 0010 0011,a0@(comm)    | mul fp0 -> fp2
  169.     movew    #0x0123,a0@(comm)        | mul fp0 -> fp2
  170.     .word    0x4a68,0xfff0,0x6bfa        | test
  171.  
  172. |    movew    #%0111 0100 1000 0000,a0@(comm)    | fetch fp1
  173.     movew    #0x7480,a0@(comm)        |
  174.     .long    0x0c688900, 0xfff067f8
  175.     movel    a0@(zahl),a1@
  176.     movel    a0@(zahl),a1@(4)
  177.  
  178. |    movew    #%0111 0101 0000 0000,a0@(comm)    | fetch fp2
  179.     movew    #0x7500,a0@(comm)        | 
  180.     .long    0x0c688900, 0xfff067f8
  181.     movel    a0@(zahl),a1@(8)
  182.     movel    a0@(zahl),a1@(12)
  183. ");    /* end asm    */
  184. #endif    sfp004
  185.  
  186.  
  187. #if defined (__M68881__) || defined (sfp004)
  188. # ifdef ERROR_CHECK    /* no error checking for now    */
  189. __asm("
  190.     pea    _funcname
  191.     jmp    c_err_check
  192. ");    /* end asm    */
  193. # else  ERROR_CHECK
  194.  
  195. __asm("rts");
  196.  
  197. # endif ERROR_CHECK
  198. #endif defined (__M68881__) || defined (sfp004)
  199.